{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 创建包含异常点的单变量数据集\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean=(-0.16) and Standard Deviation (1.04)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n_samples = 100\n",
    "fraction_of_outliners = 0.1\n",
    "number_of_inliners = int((1 - fraction_of_outliners) * n_samples)\n",
    "number_of_outliners = n_samples - number_of_inliners\n",
    "\n",
    "normal_data = np.random.randn(number_of_inliners, 1)\n",
    "mean = np.mean(normal_data)\n",
    "std = np.std(normal_data)\n",
    "print('Mean=(%0.2f) and Standard Deviation (%0.2f)' % (mean, std))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size of input data=(100,1)\n"
     ]
    }
   ],
   "source": [
    "outlier_data = np.random.uniform(low=-9, high=9, size=(number_of_outliners, 1))\n",
    "total_data = np.r_[normal_data, outlier_data]\n",
    "print('Size of input data=(%d,%d)' % (total_data.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x16bdcfc96d8>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGpFJREFUeJzt3XuwLWdd5vHvcxICbEgMOTmIJtlnE0WHDINANhhUxIFohYvEshh1ZoNxsGaX0dJAOYPEM1pg1RmUoVC8FNY2QTHZghqjRKzB4Zax0CHMPlxCMFwinnMSLmYngYAcJMTzmz+6d2VlsS7dq+/vej5Vq/Ze97dXdz/99vu+3a2IwMzM0rGv6wKYmVm9HOxmZolxsJuZJcbBbmaWGAe7mVliHOxmZolxsJuZJcbBbr0j6aiki1v4nldKurZAWb4i6UuSviDp7yT9lKRC646kNUkh6dR6Sm02n4PdbL4fjIjTgYPArwK/AFzdbZHMpnOwW69J+glJ75X0Wkmfl/SPkp4z8vyNkl4t6f2S7pX0Vkln5c99n6Q7xj7vqKSLJV0C/CLwo5L+WdKH55UlIu6NiBuAHwUuk/SE/DOfJ+mDkr4o6XZJrxx529/kf7+Qf8/TJX2LpHdLulvSXZK2JZ1Z6YcyG+FgtyH4TuDjwNnAa4CrJWnk+R8HXgJ8M3A/8JvzPjAi3g78D+CPI+KREfEdRQsTEe8H7gCekT/05bwMZwLPAy6X9EP5c9+b/z0z/57/Cwh4dV7exwPnAa8s+v1m8zjYbQiORcTvRcS/Am8Cvgn4xpHnr4mIWyLiy8AvAT8i6ZSGy/QZ4CyAiLgxIj4SEScj4mbgzcAzp70xIm6LiHdExFcjYhd43azXm5XlDh0bgs/t/RMRJ/LK+iNHnr995P9jwEPIavdNOge4B0DSd5K1vT8BOA14KPCn094o6dFkexXPAE4nq2B9vuHy2hJxjd1ScN7I/6vA14C7yJpIVvaeyGvxB0Zeu9CpTSU9lSzY35s/9EfADcB5EfENwO+SNbdM+45X548/MSLOAF408nqzyhzsloIXSbpA0grwK8B1ebPNJ4CH5Z2bDwH+O1ltes8/AWslhi6eIen5wFuAayPiI/lTpwP3RMS/SHoa8J9G3rYLnATOH3nsdOCfyTpUzwH+W9kJNpvFwW4puAb4A7Imm4cBPwfZKBbgp4GrgE+T1eBHR8nsNZfcLekDMz7/LyV9iazJ5xBZm/h/Hnn+p4FfyV/zy8Cf7D0RESeAw8Df5uPgLwJeBTwFuBf4K+D68pNsNp18oQ0bMkk3ktWer+q6LGZ94Rq7mVliHOxmZolxU4yZWWJcYzczS0wnByidffbZsba21sVXm5kN1pEjR+6KiAPzXldLsOcnMLqK7Mi7AF6SnxNjorW1NXZ2dur4ajOzpSHpWJHX1VVjfz3w9oh4oaTTGDnaz8zM2lU52CWdQXYGu58AiIj7gPuqfq6ZmS2mjs7T88kOm/79/JzUV0l6RA2fa2ZmC6gj2E8lOzz6DRHxZLLDtl8x/iJJm5J2JO3s7u7W8LVmZjZJHcF+B3BHRNyU37+OLOgfJCK2ImI9ItYPHJjbqWtmZguqHOwR8Tngdknfnj/0bODvq36umVnXtrdhbQ327cv+bm93XaJi6hoV87PAdj4i5lM8+Mx3ZmaDs70Nm5tw4kR2/9ix7D7AxkZ35Sqik1MKrK+vh8exm1mfra1lYT7u4EE4erTt0mQkHYmI9Xmv8ykFzMwmOH683ON94mA3M5tgdbXc433iYDczm+DwYVgZO4Z+ZSV7vO8c7GZmE2xswNZW1qYuZX+3tvrfcQodnd3RzGwINjaGEeTjXGM3M0uMg93MLDEOdjOzxDjYzcwS42A3M0uMg93MbMRQT/w1ysMdzcxyQz7x1yjX2M3McocOPRDqe06cyB4fEge7mVluyCf+GuVgNzPLDfnEX6Mc7GZmuSGf+GuUg93MLDfkE3+N8qgYM7MRQz3x1yjX2M3MEuNgNzNLTG3BLukUSR+U9La6PtPMzMqrs8Z+BXBrjZ9nZmYLqCXYJZ0LPA+4qo7PMzOzxdVVY/8N4OXAyWkvkLQpaUfSzu7ubk1fa2Zm4yoHu6TnA3dGxJFZr4uIrYhYj4j1AwcOVP1aMzOboo4a+3cDL5B0FHgL8CxJ19bwuWZmtoDKwR4RV0bEuRGxBvwY8O6IeFHlkpmZ2UI8jt3MrAVtXsCj1lMKRMSNwI11fqaZ2dC1fQEP19jNzBrW9gU8HOxmZg1r+wIeDnYzs4a1fQEPB7uZWcPavoCHg93MrGFtX8DDF9owM2tBmxfwcI3dzCwxDnYzs8Q42K20No+gM7Py3MZupbR9BJ2Zlecau5XS9hF0Zlaeg91KafsIOjMrz8FupbR9BJ2Zledgt1LaPoLO+skd6P3mYLdS2j6CzvpnrwP92DGIeKAD3eHeH4qI1r90fX09dnZ2Wv9eM6tubS0L83EHD8LRo22XZrlIOhIR6/Ne5xq7mZXiDvT+c7CbWSnuQO8/B7uZleIO9P6rHOySzpP0Hkm3SvqopCvqKJiZ9ZM70B+sjyOE6qix3w/8fEQ8HrgI+BlJF9TwuWaF9HHFSt3GRtZRevJk9neZQ72PI4QqB3tEfDYiPpD//yXgVuCcqp9rVkRfVyxbDn09xUatbeyS1oAnAzfV+blm0/R1xbLl0NcRQrUFu6RHAn8GvDQivjjh+U1JO5J2dnd36/paW3J9XbFsOfR1hFAtwS7pIWShvh0R1096TURsRcR6RKwfOHCgjq816+2KZcuhryOE6hgVI+Bq4NaIeF31Ig2HO+26V2bF8vyyuvV2hFBEVLoB3wMEcDPwofz23FnvufDCC2Porr02YmUlIuuyy25S9vfgwez5vrr22qyMUv/LWkSR6Zk0v1ZWhj/ttlyAnSiQyz5XzIKmnS9jz8pKT7bcY8avgAT9LWudfH4TS0HRc8U42Be0b19W75ulj6GxrAE3bX5J2VhssyHwScAaVqRzro8jM5Z1FIk7WR/M/Q3N6Mvv6mBf0KROu3F9DI1lDbi+jl7ogg/qakavftciDfF131LoPI14oNNutOO07x1zy9yJmFqn8aL2ltnx28GDXZds2Nr4XXHnabu2t7OjHY8fz2q/hw/3tzNySGW1+rm/oRlt/K5uY2/ZkE6KNKSyWv26bo7rSzt03br+XUc52K3XUg2BLrXV3zBp3vWqHbpmverHKdJeU/ctlTZ2a9Yy9wc0ren+hmnzbv/+tNv3m/5dSbmNfZE2YrcrD8+yjrlPwbwD+Ma5fb+YZNvYZ+3KTdttT3n3L2XLOuY+BWXnUerDbds2uGCfdv7tK66YHt59PGe3247n61NnVMqaWBanzaP9+3vUDp2yIu01dd+qtLGPjxefd9tr75r0nLRwMSqZ13bs8dYZt7E3r6nfeNbnll2+vT48gIJt7IML9mkHAUy77S0MfeqwmVWeZQ2zaSvvIiu1g6C4JteNomfdnPWaZV0fpkk22CfN6Hk19i4XjkkL7qw9iL5thNpQdP749Lz163Jvtsi8Wsb1YZZkgz3iwYfyz7p13byxyJCvvjUbtaHIyls0sGctF32ovfdtb6LL4Czy3cu4PsySdLDvmdXevn9/dutiBZq34dm/f3pILWMNpcjKW/R3mdcH02XtvY97E12Wqc75viyWItinzfRZwTlLHbWpIk1F0uw25b6t/E2rUnMbr4kX2ZPrY99Kl7rai6hzT21ZLEWw13l0W10LUB3B0rfd9aIWLXeVttbx1xfdsHbBzQoPVmffyrJYimCPKN85OU1dtak+NwU0oa5TF4/Ox0nNaEUCe29ezWsKc419urZD1EMfy1maYJ9kkRWortpU3zvv6lQmbKe9f3wlLTL+uWhNvG+78X0rz7i+H1/R99+vDa0GO3AJ8HHgNuAV817fdLAvsgDUVZtapoWvSLPTtA1jlWa0MvOqjrHURV9TRNfhOEvfj6+oo9LU59+/iNaCHTgF+AfgfOA04MPABbPes0iwN7XLVldTQpWyDlWRo4CnbRiLbBSmbSDqDJkin9WHUBstS1PLVt+Pr6jazNmn+bioNoP96cBfj9y/Erhy1nvKBntTM2TS5+4tPCkHcl2KdmhOssipIUbVFXBFAqsvzWtNB9Os36JqU2Ud86vqwIQ+bJyqajPYXwhcNXL/xcBvT3jdJrAD7KyurpaamKZmSAozepam9xyqbBjrHqq6qCLDKPvSId708jprw1Hlu+sa/VJ1xFMKo5LaDPb/MCHYf2vWe8rW2JuaIUOY0Yu2/7a129nEEMc2m7KK7HVMa/dvuzLQxvI67bcvOr8mjWiqc7x6lRFPKVTkkmqKWdYae5X23yFcqaYPfRFFaoGT9iLKhGubzUbzprVKOYpWIMaXxSK/WdlpW6Ti4jb2csF+KvAp4LEjnaf/dtZ7+tzG3qcZXbX9t+na3SL6EObjigyjXLSm2HZHbxPvnaXI8nfKKfN/s0X2RhZZlvq4/JXR9nDH5wKfyEfHHJr3+jZGxdT5uV0tDEUW9qqdkG1oYuRRE+psMijzmWUsuiw2tXdadPnzWRzrsdQHKNWpy1p9lRp7nZ2QVTZsRZo6+rLy1tXJN2rexrmtSkNT7fNFR6os0jHap41+XzjYa9JlTaLqGOs6QqPqCldkxe+6eWhU1Y3YeCfivLBrK8yaWo6LtLGX6VAfcjNJGxzsNamzptNUm2BdYTTpvVUDocpBTENSZM9kPOzmtdnXGWxNbkTmjYqx+jjYZygThHXVdPq4q1mkTFU3bEWGE6aw8hftxB5d3toeH+8a8fA52KcoG7B1BXLXnUOTVuoqbfhlhtgtw9G9RfZMxjeGRduny3KAp8vBPsUiQVXHitJlJ9q0jVORAFp0w7Zsu+eLhHSR5puyTX593DO0+jjYp+jqaNNZG5SmV8Zp311kfHFE+Y3OMobLop2Ii46Pn6brPcO6eK9jMgf7FEUX/GkL1qIL3Kywa3plnNVMMK2ZpEotO5VwKavKXkpdG8MhnCZjnmWsGBTlYJ+iyhDCyy+vtsBN2yg0vTLO21uYdvDQotOZQrh0oY5aagob1RSmoSkO9hkWHeJXtOmirKYX5CIbszo78rxidieFzmpXDKZzsFdQ9jD9qgtcG7ue8zZmi4zq6HJ6+qKPbcGz9sKGMB9cMZjOwV5B2zX2iO4Dou6hd11PTxu63oA1fXBZV+r+XVNaFh3sFTTVxt5ndR4aviyaDM55HbFtHFzWpbrCuOuNb90c7BXVPSqmrfLV9ZnLMPa8qqaCs8hGtu6Lfqcqtd/AwZ6w1GohXejzCJSipyeYt1FZtuVk0jwd8l7LJA72hKVWC2lbXYHXVHCW7byftQxU3YD1ZQ91nmnzYghXEivDwZ6w1Gohbatzw9hE8BWpsbdx0e8h1fin/WZtXxy9aUWDXdlr27W+vh47Ozutf28q1tbg2LGvf/zgQTh6tO3SDM++fdkqPk6CkyfbL8+47W3Y3IQTJyY/v7ICW1vZ/4cOwfHjcNZZ2f177oHVVTh8GDY2qpVjSMvZtHkKsH9/9rfO36Yrko5ExPq81+1rozBWr8OHs5V71MpK9rjNt7pa7vG2bWxkwX3wYLax2b8/u0nZY1tb2Ws2NrKAveYa+MpX4O67s3A7dizbMGxvVyvH8ePlHu/SrHl3993Z73PNNdnvNdRQL6VItb7um5tiqhtK22cfDamJoYi2O3H72D5d5EyZfSx3WRRsiqlUY5f0PyV9TNLNkv5c0pk1bW9sjr3a2smTS1QLqcl4jXi0FjxETdWsh7RnODpPp+njnkZTqjbFvAN4QkQ8EfgEcGX1Ipk1L6UNY1NNS0PbAO7N02nh3pemtjZUCvaI+N8RcX9+933AudWLZGZlNFmzHuIGcEh7Gk2ps/P0JcD/qvHzbMz2djZSYd++7G/VzjFLw9Bq1k3z78H84Y6S3gk8ZsJThyLirflrDgHrwA/HlA+UtAlsAqyurl54bNI4Kptq0hC4vWFvy7TAmi2zosMdK49jl3QZ8FPAsyNiysjbB/M49vKGNKbYzJpRNNhPrfgllwC/ADyzaKjbYoY0ptjMulW1jf23gdOBd0j6kKTfraFMNkHfD6oxs/6oVGOPiG+tqyA22+HDk9vYl6mn38yK8SkFBsI9/WZWlIN9QIqMKfaQSDNzsCdkb0jksWP1ngzKLEUpV4Ic7Ak5dOjrT/V64kT2uJk9IPVKkIM9IR4SWa+Ua3TLLvVKkIM9IR4SWZ/Ua3TLLvVKkIM9IT75UX1Sr9Etu9QrQQ72hHhIZH1Sr9Etu9QrQQ72xAzxNKt9lHqNbtlNuvzgwx8OL35xGv0pDnazCVKv0Vnz14ztkoPdbAI3ay2PFPtTKp+2dxE+ba+Z9cW+fVlNfZyUNWn2SdHT9rrGbmZLLcX+FAe7mS21FPtTHOxmttRS7E+pdD52M7MUbGwMO8jHucZuZpYYB7uZWWIc7GZmiXGwm5klppZgl/RfJYWks+v4PDMzW1zlYJd0HvD9gM97Z2bWA3XU2H8deDnQ/rkJzMzs61QKdkkvAD4dER8u8NpNSTuSdnZ3d6t8rZmZzTD3ACVJ7wQeM+GpQ8AvAj9Q5IsiYgvYguwkYCXKaGZmJcwN9oi4eNLjkv4d8Fjgw5IAzgU+IOlpEfG5WktpZmaFLXxKgYj4CPDovfuSjgLrEXFXDeUyM7MFeRy7mVlJ29vZJfT27evnpfRqOwlYRKzV9VlmZn21vZ1dOm/vqkt7l9KD/pxIzDV2M7MShnApPQe7mVkJx6ccijnt8S442M3MShjCpfQc7GZmJQzhUnoOdjOzEoZwKT1fGs/MrKS+X0rPNXYzs8Q42M3MEuNgNzNLjIPdzCwxDnYzs8Q42M3MEuNgNzNLjIPdzCwxDnYzs8Q42M3MEuNgNzNLjIPdzCwxDnYzs8RUDnZJPyvp45I+Kuk1dRTKzMwWV+m0vZL+PXAp8MSI+KqkR9dTLDMzW1TVGvvlwK9GxFcBIuLO6kUyM7Mqqgb7twHPkHSTpP8j6anTXihpU9KOpJ3d3d2KX2tmZtPMbYqR9E7gMROeOpS//1HARcBTgT+RdH5ExPiLI2IL2AJYX1//uufNzKwec4M9Ii6e9pyky4Hr8yB/v6STwNmAq+RmZh2p2hTzF8CzACR9G3AacFfVQpmZ2eKqXsz6jcAbJd0C3AdcNqkZxszM2lMp2CPiPuBFNZXFzMxq4CNPzcwS42A3M0uMg93MLDEOdjOzxDjYzcwS42A3M0uMg93MrAHb27C2Bvv2ZX+3t9v77qoHKJmZ2ZjtbdjchBMnsvvHjmX3ATY2mv9+19jNzGp26NADob7nxIns8TY42M3Manb8eLnH6+ZgNzOr2epqucfr5mA3M6vZ4cOwsvLgx1ZWssfb4GA3M6vZxgZsbcHBgyBlf7e22uk4BY+KMTNrxMZGe0E+zjV2M7PEONjNzBLjYDczS4yD3cwsMQ52M7PEONjNzBJTKdglPUnS+yR9SNKOpKfVVTAzM1tM1Rr7a4BXRcSTgF/O75uZWYeqBnsAZ+T/fwPwmYqfZ2ZmFVU98vSlwF9Lei3ZRuK7pr1Q0iawCbDa1plwzMyW0Nwau6R3Srplwu1S4HLgZRFxHvAy4OppnxMRWxGxHhHrBw4cqG8KzMxa0OUVkcpSRCz+Zule4MyICEkC7o2IM+a9b319PXZ2dhb+XjOzNo1fEQmyszW2eWIvAElHImJ93uuqtrF/Bnhm/v+zgE9W/Dwzs97p+opIZVVtY/8vwOslnQr8C3kbuplZSrq+IlJZlYI9It4LXFhTWczMeml1Nbsg9aTH+8hHnpqZzdH1FZHKcrCbmc3R9RWRyvIVlMzMCujyikhlucZuZpYYB7uZWWIc7GZmiXGwm5klxsFuZpaYSueKWfhLpV1gwnD/Qs4G7qqxOEOxjNO9jNMMyzndyzjNUH66D0bE3LModhLsVUjaKXISnNQs43Qv4zTDck73Mk4zNDfdbooxM0uMg93MLDFDDPatrgvQkWWc7mWcZljO6V7GaYaGpntwbexmZjbbEGvsZmY2g4PdzCwxgwp2SZdI+rik2yS9ouvyNEHSeZLeI+lWSR+VdEX++FmS3iHpk/nfR3Vd1rpJOkXSByW9Lb//WEk35dP8x5JO67qMdZN0pqTrJH0sn+dPT31eS3pZvmzfIunNkh6W4ryW9EZJd0q6ZeSxifNWmd/Ms+1mSU+p8t2DCXZJpwC/AzwHuAD4j5Iu6LZUjbgf+PmIeDxwEfAz+XS+AnhXRDwOeFd+PzVXALeO3P814Nfzaf488JOdlKpZrwfeHhH/BvgOsulPdl5LOgf4OWA9Ip4AnAL8GGnO6z8ALhl7bNq8fQ7wuPy2CbyhyhcPJtiBpwG3RcSnIuI+4C3ApR2XqXYR8dmI+ED+/5fIVvRzyKb1TfnL3gT8UDclbIakc4HnAVfl90V2gfTr8pekOM1nAN8LXA0QEfdFxBdIfF6TXQfi4fm1kleAz5LgvI6IvwHuGXt42ry9FPjDyLwPOFPSNy363UMK9nOA20fu35E/lixJa8CTgZuAb4yIz0IW/sCjuytZI34DeDlwMr+/H/hCRNyf309xfp8P7AK/nzdBXSXpESQ8ryPi08BrgeNkgX4vcIT05/WeafO21nwbUrBrwmPJjtWU9Ejgz4CXRsQXuy5PkyQ9H7gzIo6MPjzhpanN71OBpwBviIgnA18moWaXSfI25UuBxwLfDDyCrBliXGrzep5al/chBfsdwHkj988FPtNRWRol6SFkob4dEdfnD//T3q5Z/vfOrsrXgO8GXiDpKFkT27PIavBn5rvrkOb8vgO4IyJuyu9fRxb0Kc/ri4F/jIjdiPgacD3wXaQ/r/dMm7e15tuQgv3/AY/Le89PI+twuaHjMtUub1u+Grg1Il438tQNwGX5/5cBb227bE2JiCsj4tyIWCObr++OiA3gPcAL85clNc0AEfE54HZJ354/9Gzg70l4XpM1wVwkaSVf1vemOel5PWLavL0B+PF8dMxFwL17TTYLiYjB3IDnAp8A/gE41HV5GprG7yHbBbsZ+FB+ey5Zm/O7gE/mf8/quqwNTf/3AW/L/z8feD9wG/CnwEO7Ll8D0/skYCef338BPCr1eQ28CvgYcAtwDfDQFOc18GayfoSvkdXIf3LavCVrivmdPNs+QjZqaOHv9ikFzMwSM6SmGDMzK8DBbmaWGAe7mVliHOxmZolxsJuZJcbBbmaWGAe7mVli/j/AKjQQq6hf1AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x16bdbfa6eb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.cla()\n",
    "plt.figure(1)\n",
    "plt.title('Input Data')\n",
    "plt.scatter(range(len(total_data)), total_data, c='b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-11T16:06:05.105754Z",
     "start_time": "2020-04-11T16:06:05.100497Z"
    }
   },
   "source": [
    "# 绝对中位差检查异常点"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Median absolute Deviation=1.04\n",
      "Outlier 4.69\n",
      "Outlier 5.59\n",
      "Outlier -8.27\n",
      "Outlier -4.47\n",
      "Outlier -7.04\n",
      "Outlier 3.50\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHVJJREFUeJzt3X2UJXV95/H3Z3ok2iiCwygPw3RDonHRuCKtohJjFHOQqHg8uuvuaEg2u52dyTHoiWGFybp6TiYPrmeVRJlsR00M06smSIxxNxofk8xuxPQIKjg+Z2boADKAIjKEh5nv/lHVcKe5t2/dW1W36v7u53XOPd23bt2qXz1961ff36/qKiIwM7N0rGu6AGZmVi0HdjOzxDiwm5klxoHdzCwxDuxmZolxYDczS4wDuzVO0qykkLQ+f//Xki5quly9tL18vUj6vKT/2HQ5rH4O7DYwSb8o6auSDkm6RdJOSccP8P19ks7r9XlEvCQiPlBNaavX9vKZObDbQCT9OvB7wG8AjwXOAWaAT0k6puGyrW9y/mZt4cBuhUk6Dngb8PqI+ERE3B8R+4B/QxbcX5uP9yeSfqvjey+QtJz/fyWwGfgrST+SdEmX+RyVMpD0HyTtlfR9SZ+UNNPxWUj6VUnfAr6lzDsl3SrpTklfkfTUHstz1JWDpLdK2pX//0hJuyTdLukHkv5R0hNWly+/etkt6R15+f5J0ks6pnm6pL+TdJekT0t6z8o8upTnBZKWJV2Sl/9mSa+QdIGkb0q6Q9JlHeM/S9I/5OW7WdK7O0+ukl4s6ev5eng3oG7ztfQ4sNsgngs8Eri6c2BE/Aj4a+DF/SYQEa8DDgAvi4hHR8Tb1xpf0iuAy4BXAhuBvwc+uGq0VwDPBs4Efg54PvAk4Hjg3wK39ytXFxeRXZGcBmwA/jNwT49xnw18AzgReDvwPkkrQfR/AV/Mp/FW4HV95nsS2To+FXgL8EdkJ8yzgZ8G3iLpjHzcw8Ab8/k+B3gRsA1A0onAR4DfzD//DvC8Igtu48+B3QZxInBbRDzQ5bOb88+r9ivA70TE3ny+vw08vbPWnn9+R0TcA9wPPAZ4MqD8ezcPMd/7yYLxT0TE4YjYExE/7DHu/oj4o4g4DHwAOBl4gqTNwDOBt0TEfRGxG/hYgfnuiIj7gQ+RrdPLI+KuiLgBuAF4GkBepi9ExAP5ldP/BH4mn84FwNci4qp8Wu8CbhliPdgYcmC3QdwGnNgjl31y/nnVZoDL83TDD4A7yFIKp3aMc+PKPxHxWeDdwHuA70layFNIg7oS+CTwIUk3SXq7pEf0GPfBgBkRh/J/Hw2cAtzRMeyosvZwe36CgIeuEL7X8fk9+bSR9CRJH88bsH9IdtJbObmewtHrJQrM2xLhwG6D+AfgXrK0yIMkHQu8BPhMPuhuYLpjlJNWTWeQR4reCPxKRBzf8XpURPy/XtOLiN+PiLOBp5ClZH6jx7R7ljNvP3hbRJxJloJ6KfALA5QbsquYx0nqnMdpA05jLTuBrwNPjIjjyFJWKymgmzvnlaeGqpy3tZgDuxUWEXeSNZ7+gaTzJT1C0izw58AyWS0X4DrgAkmPk3QS8IZVk/oecAbF/CFwqaSnAEh6rKRX9xpZ0jMlPTuvXd8N/AtZLrqb64DX5MsxB7yqYzo/K+mnJE0BPyRLkfSaTlcRsR9YAt4q6RhJzwFeNsg0+nhMXrYfSXoysLXjs/8NPEXSK/MrrF/j4SdYS5QDuw0kb+y8DHgHWVC5hqxW/aKIuDcf7Urgy8A+4G+AD6+azO8Av5mnV97UZ35/Qda98kN5uuF6squDXo4ja3D8PrCfrOH0HT3G/a/Aj+fjvo2soXPFScBV+TLuBf4W6NqbpY8tZA2btwO/RbYu7l3zG8W9Cfj3wF1ky/zgeo6I24BXA7+bz/uJwP+taL7WcvIPbZiNjqQPA1+PiP/WdFksXa6xm9UoTw39uKR1ks4HLgQ+2nS5LG2+U8+sXieR9fvfQNYOsTUirm22SJY6p2LMzBLjVIyZWWIaScWceOKJMTs728SszczG1p49e26LiI39xmsksM/OzrK0tNTErM3Mxpak/UXGcyrGzCwxDuxmZompJLBLOl7SVfmzn/fmt06bmVkDqsqxXw58IiJelT/of7rfF8zMrB6lA3v+SNTnA78IEBH3AfeVna6ZmQ2nilTMGcBB4I8lXSvpvfljXI8iaV7SkqSlgwcPVjBbMzPrporAvh54BrAzIs4ie1Tqm1ePFBELETEXEXMbN/bthmlmZkOqIrAvA8sRcU3+/iqyQG9mZg0oHdgj4hbgRkk/mQ96EfC1stM1M2vC7m2LLK+f5YjWsbx+lt3bFpsu0sCq6hXzemAx7xHzXeCXKpqumdnI7N62yFk75zmW7GdqNx3ezwk759kNnHvFlmYLN4BGnu44NzcXfqSAmbXN8vpZNh1++F37y1MzbHpg3+gLtIqkPREx128833lqZpY75fCBgYa3lQO7mVnupqnNAw1vKwd2M7Pcvvkd3L3qxvm7mWbf/I6GSjQcB3Yzs9y5V2zh2q0LLE/NcASxPDXDtVsXxqrhFNx4amY2Ntx4amY2oRzYzcwS48BuZpYYB3Yzs8Q4sJuZJcaB3cwmXgoP/upU1UPAzMzGUioP/urkGruZTbTZhe0PBvUVx3KI2YXtDZWoPAd2M5toqTz4q5MDu5lNtFQe/NXJgd3MJloqD/7q5MBuZhMtlQd/dfJDwMzMxoQfAmZmNqEqC+ySpiRdK+njVU3TzMwGV2WN/WJgb4XTMzOzIVQS2CVtAn4eeG8V0zMzs+FVVWN/F3AJcKTXCJLmJS1JWjp48GBFszUzs9VKB3ZJLwVujYg9a40XEQsRMRcRcxs3biw7WzMz66GKGvvzgJdL2gd8CHihpF0VTNfMzIZQOrBHxKURsSkiZoHXAJ+NiNeWLpmZmQ3F/djNzGrS1HPeK30ee0R8Hvh8ldM0MxtHTT7n3TV2M7MaNPmcdwd2M7MaNPmcdwd2M7MaNPmcdwd2M7MaNPmcdwd2M7MaNPmcdz+P3cxsTPh57GZmE8qB3cwsMQ7sNpDFRZidhXXrsr+Lo7mRzswGUOmdp5a2xUWYn4dD+T0X+/dn7wG2jO/v/polxzV2K2z79oeC+opDh7LhZtYeDuxW2IEeN8z1Gm5mzXBgt8I297hhrtdwM2uGA7sVtmMHTB99Ix3T09lwS58bzseHA7sVtmULLCzAzAxI2d+FBTecToKVhvP9+yHioYZzB/d28p2nZtbX7GwWzFebmYF9+0ZdmsnlO0/NrDJuOB8vDuxm1pcbzseLA7uZ9eWG86O1vSG5dGCXdJqkz0naK+kGSRdXUTCzItp+gKXCDecPGYeG5NKNp5JOBk6OiC9JegywB3hFRHyt13fceGpVWP2IA8hqkZMacGw0mmxIHlnjaUTcHBFfyv+/C9gLnFp2umb9+BEH1oRxaEiuNMcuaRY4C7imy2fzkpYkLR08eLDK2dqEGocDzNIzDg3JlQV2SY8GPgK8ISJ+uPrziFiIiLmImNu4cWNVs7UJNg4HmKVnHBqSKwnskh5BFtQXI+LqKqY5Dtxw16yiB5i3k1VpLBqSI6LUCxDwp8C7in7n7LPPjnG3a1fE9HRE1i6evaans+FttmtXxMxMhJT9bXt5++m3POO6ncy6AZaiQIytolfMucDfA18FjuSDL4uI/9PrOyn0iunVMg7ZGXzHjpadwZnMXiS+Fd5SUrRXjJ8VM6R167L6Xy9tDJiTGOR6bScJjhx5+HCzNvOzYmrWr4Gujd3uJrEXiRtYj+b2huq1cZ06sA+pW8Pdam0LmJMY5MahB8OojMMdk+OmrevUgX1InS3jvbQtYE5ikBuLHgwj4hu6qtfWdeocewXGqVFycTHb6Q4cyE48bWzktXq4vaF6o16nzrGP0DjVCrdsyRpKjxzJ/raxjFaPplNxbcxFl9X0Ou3Fgb0iDpjVSzEQNGlUqbhu262tueiyWpveLNLZverXqG9QSu2mnEngG4vqUfex0Gu7bdhw9LCV18xMtfNvwijjC6O6QWkYZXPsg+SJxyn/bQ+ZxD73KVjrxr1unN8fTLI59rUu6bpdAra11drWNol97lMw6PZpOhedqrEL7L0C9cUXdw/4vWoPTQcI54/X1tZGqZTUsQ/22j4bNrQ0F52qIvmaql9lcuxS91xdr9fUVPtye2vlj90ekBkmx+51V1xdbRhV7tveng9HwRz72AX2mZnBAvvKjtVUI1y3nbPXMmzYMJkNhr0O4M7hGzZkLz/FsRq99sEqKjxFArKfyjmcZAN7tw2+1qszmI76zN9r5xz0xJRCz4FeihzARcZZ64Tv2t7D9bryleqfd5ntmfKxUESygT1i7Vpvm87wvcrYKz3U6zWKg60pRQ7gIuP0S9E1vS9EtCu10GTgLLM9Uz4Wikg6sK9Y62Dud+neTVUHXpETz6T19e2lyAFcZJwiJ/q2tqtMWnnKbM+Uj4UiJiKwV5mrrmpHL5Iq6pUeatvBP4hhT4plauyd67LIem+yttfGQNXUFUSRdTHOx0KdJiKwV3mXW1UHXpGaemo9O8ochMPm2LuN3+9Kqckg2vbUwij3u6L7yzgeC3WbiMAe0X3jD3MQVXXgrZUeSm3nrCqQFun9Msi82ljba2ONfUUT3W/d9XE4Iw3swPnAN4BvA2/uN37dz4oZ5iCqu8behgO4SmVTH8OmooqegNsWCNp4slnR9u63Zddd2/aFMkYW2IEp4DvAGcAxwJeBM9f6Tt2BfdibW+rKsbflAK5SmcbKMim0Kk+cVfS3rnp+TRj0pr9RV1LKdGVN7XgcZWB/DvDJjveXApeu9Z1hAntdl26D3ARTV1nHUZnuhYPeZNZZGx/lCbhNQaHOfarM9hhFuevY18b1CnqUgf1VwHs73r8OePda3xk0sNd1gLXpwB03ZWpRZWuIVQSLMr1xpqZGe9Kuez+tshNCHeUuc3XY9kbrQY0ysL+6S2D/gy7jzQNLwNLmzZsHWpi6zrqpnc071X3VUOagbUNOt0gjd5ET0CgqAqPYT8t2v+32/aLl7revlmnPSe0YTyoVU9dZt+1n82FzwKO6Ehn25NFEL4zVinRL7VVjHSRI1JmKGNXt/73K3xm8V5dxrUA8TGpt2B5YqV2VjzKwrwe+C5ze0Xj6lLW+M2hgn8Qae5kc8Djcwdp0O0SRWmC3K4hBAmzdqYhhupNWta6LrL8iT1YddNmG7RiRSpvXqLs7XgB8M+8ds73f+OOWY29ixyiTAx40ANWpjsbpqsu21vrqLP+gj4CuquJQZv+v69gpuu9V1X119TKlEqgHldwNSnVtzGHye6O4lCvzDJVer1HX2PvV6tpySTxILniQfaHKFMqw+39dV6VF9r3OXHuvcrf5qrmNkgvsTWlqxytTY6+6EbLqoNK2A7hsI2Eva23DUdU668rPF2mjaPKKIlUO7BVpquGqbD/rYRtehylHL0VqdePUUF30+yspp5XlW73utm4dXTCrq2LSbb9YWdZB198kp1YG5cBekSoPjEF34KqCc69pFwkuZZZ/XGrsZfVLOa0OeP16d1QZ2OqsETsgj54Dew/DBNdxf9RAmT7GZa5YxiXHXtagJ7C6fhSk177tAJwOB/Yuhg2uVRwYTeVbey1z0dTIMDX2Xj1h2tYrpiqDppzquJJxrnoyOLB30WQL/FoHf50HZK9lLtp1b9CAMYkBZtBAXaQP+KBtD6n0LvHVxdoc2LuooiG06h4ig/aNHtQgJ5TOXPDWrcP1P08lwAximJRTkVz7INp+F3URk1gpGJQDexeDBJ1uAbyOG0WqqrENs8ydwaWqvG8KAWYYw96IVVUwS+GEmsIy1M2BvYtB7jSt62l3wzZi1rnMRe8iLNMTxgdnb1WkH6rsftiUSa0UDMKBvYciB1HRQFfFjjeKy89+y1z0DtZhe8KkejndtnzwWldg47ANXCnoz4G9hFHfqt90gKiyxh7R/PKMQtu6r3Ya1wBZ5TpNdR90YC+h14HRlt+ArFqRXhopLGeV6ryjc61cfZHgN84pjbrSUqnsvw7sJay1Y7ShJlBHGVZPs7NXTEo1nqrWXR3Bs0jvmlH/Luw4Snn5HdhLakMA7ybl2kjdqlx3dQSPQdt2ep1QJmkf6XacjvMVSz8O7IlKuTZStyrXXR3Bc9C2nbWWoUzFpK2VmtXq6r3WZg7siUq5NlK3qtdd1QGwSI297naecartr7W+xrFXUBFFA/s6bKxs3jzYcHtI1etuyxbYtw+OHMn+btkybMkyO3bA9HTvz6en4fLLYWEBZmZAgg0b4FGPgte9DmZnYXGxXBm2b4dDh44eduhQNrxtDhzo/VlEtn4gW1cLC+W3zzhxYB8z3Q7+6elsuK2t7etuy5aHB+0NG7L/O4PTygnlyivhnnvg9tuzQLZ/P8zPlwvuvYLlWkG0Kf1OyBHZeqvipDt2ilTrq345FVPOuORA2yildTfKBtw25qfreJha21EwFaNs3OFI+u/Ay4D7yH7I+pci4gf9vjc3NxdLS0tDz9fMYN26LHytJmXpoWEsLma1/s50zPR0e1MZi4tZmmj//u6fr9TYUyFpT0TM9RuvbCrmU8BTI+JpwDeBS0tOz8wKqqO9ZXU6qO356ZW01K5d7U6zjVqpwB4RfxMRD+RvvwBsKl8k62ZxMWscW7eumkYyG391tRlU3Sg8CuN2QqpbqVTMUROS/gr4cETs6vH5PDAPsHnz5rP397p2socZt8tjG52VVMSBA1lNfccO7xMpK5qK6RvYJX0aOKnLR9sj4i/zcbYDc8Aro8CZwjn2wczOds8hppY/NLO1FQ3s6/uNEBHn9ZnRRcBLgRcVCeo2uHHqgmZmzSuVY5d0PvBfgJdHxKF+49twfFOSmQ2ibK+YdwOPAT4l6TpJf1hBmWyVtt9YY2bt0jcVs5aI+ImqCmK9rTSGuZHMzIooFdhtdFZuJTcz68fPikmI+7qbrW1SjhHX2BOxuq/7ygOhwDV9M5isY8Q19kSM0+NW225SanWTZpKOEdfYE+G+7tWYpFrdpJmkY8Q19kS4r3s1JqlWN2km6RhxYE+E+7pXY5JqdZNmko4RB/ZE+Ol21ZikWt2k6fYLVVX+rGCbOLAnZBwft9o2k1Srm0R1/qxgmziwm3Xwlc9kSL0tpbLnsQ/Cj+01sybV8bOCozCqn8YzMxs7qbelOLCb2cRJvS3Fgd3MJk7qbSm+89TMJlLKT0x1jd3MLDEO7GZmiXFgNzNLjAO7mVliKgnskt4kKSSdWMX0zMxseKUDu6TTgBcDfv6dmVkLVFFjfydwCTD6ZxOYmdnDlArskl4O/HNEfLnAuPOSliQtHTx4sMxszcxsDX1vUJL0aeCkLh9tBy4Dfq7IjCJiAViA7CFgA5TRzMwG0DewR8R53YZL+ingdODLkgA2AV+S9KyIuKXSUpqZWWFDp2Ii4qsR8fiImI2IWWAZeIaDupmlaPe2RZbXz3JE61heP8vube39VQ4/K8bMrI/d2xY5a+c8x5L9Osemw/s5Yec8u4Fzr2jfA2cqu0Epr7nfVtX0zMzaYnZh+4NBfcWxHGJ2oZ0/ueQ7T83M+jjlcPfbdHoNb5oDu5lZHzdNdf9ppV7Dm+bAbmbWx775HdzN0T+5dDfT7Jtv508uObCbmfVx7hVbuHbrAstTMxxBLE/NcO3WhVY2nAIouv1Ud83m5uZiaWlp5PM1MxtnkvZExFy/8VxjNzNLjAO7mVliHNjNzBLjwG5mlhgHdjOzxDiwm5klxoHdzCwxDuxmZolxYDczS4wDu5lZYhzYzcwS48BuZpYYB3Yzs8SUDuySXi/pG5JukPT2KgplZmbDK/Vj1pJ+FrgQeFpE3Cvp8dUUy8zMhlW2xr4V+N2IuBcgIm4tXyQzMyujbGB/EvDTkq6R9LeSnllFoczMbHh9UzGSPg2c1OWj7fn3TwDOAZ4J/JmkM6LLzzJJmgfmATZvbucPwJqZpaBvYI+I83p9JmkrcHUeyL8o6QhwInCwy3QWgAXIfhpv6BKbmdmayqZiPgq8EEDSk4BjgNvKFsrMzIZXqlcM8H7g/ZKuB+4DLuqWhjEzs9EpFdgj4j7gtRWVxczMKuA7T83MEuPAbmaWGAd2M7PEOLCbmSXGgd3MLDEO7GZmiXFgNzOr0O5tiyyvn+WI1rG8fpbd2xZHXoayNyiZmVlu97ZFzto5z7EcAmDT4f2csHOe3cC5V2wZWTlcYzczq8jswvYHg/qKYznE7ML2kZbDgd3MrCKnHD4w0PC6OLCbmVXkpqnujyTvNbwuDuxmZhXZN7+Du5k+atjdTLNvfsdIy+HAbmZWkXOv2MK1WxdYnprhCGJ5aoZrty6MtOEUQE08ZXdubi6WlpZGPl8zs3EmaU9EzPUbzzV2M7PEOLCbmSXGgd3MLDEO7GZmiXFgNzNLTKnALunpkr4g6TpJS5KeVVXBzMxsOGVr7G8H3hYRTwfekr83M7MGlQ3sARyX//9Y4KaS0zMzs5LKPrb3DcAnJb2D7CTx3F4jSpoH5gE2bx7tcxPMzCZJ3xq7pE9Lur7L60JgK/DGiDgNeCPwvl7TiYiFiJiLiLmNGzdWtwRmZjVpw49mDKPUIwUk3QkcHxEhScCdEXFcv+/5kQJm1narfzQDsgd6NfHslxWjeqTATcDP5P+/EPhWyemZmbVCW340Yxhlc+z/Cbhc0nrgX8hz6GZm464tP5oxjFKBPSJ2A2dXVBYzs9a4aWozmw7v7z68gfIMwneempl10ZYfzRiGA7uZWRdt+dGMYfiHNszMxoR/aMPMbEI5sJuZJcaB3cwsMQ7sZmaJcWA3M0uMA7uZWWIc2M3MEuPAbmaWmEZuUJJ0EHj4QxiKORG4rcLijAMv82TwMk+GMss8ExF9f9CikcBehqSlIndepcTLPBm8zJNhFMvsVIyZWWIc2M3MEjOOgX2h6QI0wMs8GbzMk6H2ZR67HLuZma1tHGvsZma2Bgd2M7PEjFVgl3S+pG9I+rakNzddnjpIOk3S5yTtlXSDpIvz4Y+T9ClJ38r/ntB0WaskaUrStZI+nr8/XdI1+fJ+WNIxTZexSpKOl3SVpK/n2/o5E7CN35jv09dL+qCkR6a2nSW9X9Ktkq7vGNZ1uyrz+3k8+4qkZ1RVjrEJ7JKmgPcALwHOBP6dpDObLVUtHgB+PSL+FXAO8Kv5cr4Z+ExEPBH4TP4+JRcDezve/x7wznx5vw/8ciOlqs/lwCci4snAvyZb9mS3saRTgV8D5iLiqcAU8BrS285/Apy/aliv7foS4In5ax7YWVUhxiawA88Cvh0R342I+4APARc2XKbKRcTNEfGl/P+7yA74U8mW9QP5aB8AXtFMCasnaRPw88B78/cCXghclY+S2vIeBzwfeB9ARNwXET8g4W2cWw88StJ6YBq4mcS2c0T8HXDHqsG9tuuFwJ9G5gvA8ZJOrqIc4xTYTwVu7Hi/nA9LlqRZ4CzgGuAJEXEzZMEfeHxzJavcu4BLgCP5+w3ADyLigfx9atv6DOAg8Md5+um9ko4l4W0cEf8MvAM4QBbQ7wT2kPZ2XtFru9YW08YpsKvLsGT7akp6NPAR4A0R8cOmy1MXSS8Fbo2IPZ2Du4ya0rZeDzwD2BkRZwF3k1DapZs8r3whcDpwCnAsWSpitZS2cz+17efjFNiXgdM63m8CbmqoLLWS9AiyoL4YEVfng7+3cpmW/721qfJV7HnAyyXtI0uvvZCsBn98fskO6W3rZWA5Iq7J319FFuhT3cYA5wH/FBEHI+J+4GrguaS9nVf02q61xbRxCuz/CDwxb0U/hqzh5WMNl6lyeX75fcDeiPgfHR99DLgo//8i4C9HXbY6RMSlEbEpImbJtulnI2IL8DngVfloySwvQETcAtwo6SfzQS8Cvkai2zh3ADhH0nS+j68sc7LbuUOv7fox4Bfy3jHnAHeupGxKi4ixeQEXAN8EvgNsb7o8NS3juWSXY18BrstfF5DlnT8DfCv/+7imy1rDsr8A+Hj+/xnAF4FvA38O/FjT5at4WZ8OLOXb+aPACalvY+BtwNeB64ErgR9LbTsDHyRrQ7ifrEb+y722K1kq5j15PPsqWY+hSsrhRwqYmSVmnFIxZmZWgAO7mVliHNjNzBLjwG5mlhgHdjOzxDiwm5klxoHdzCwx/x87Bmc+WUMzHAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x16bdd018438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "median = np.median(total_data)\n",
    "b = np.sqrt(2)\n",
    "mad = b * np.median(np.abs(total_data - median))\n",
    "outliers = []\n",
    "outlier_index = []\n",
    "print('Median absolute Deviation=%0.2f' % (mad))\n",
    "lower_limit = median - (3 * mad)\n",
    "upper_limit = median + (3 * mad)\n",
    "for i in range(len(total_data)):\n",
    "    if total_data[i] < lower_limit or total_data[i] > upper_limit:\n",
    "        print('Outlier %0.2f' % (total_data[i]))\n",
    "        outliers.append(total_data[i])\n",
    "        outlier_index.append(i)\n",
    "\n",
    "plt.figure(2)\n",
    "plt.title('Outliers using mad')\n",
    "plt.scatter(range(len(total_data)), total_data, c='b')\n",
    "plt.scatter(outlier_index, outliers, c='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 平均值加减3倍标准差检查异常点\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Outlier 5.59\n",
      "Outlier -8.27\n",
      "Outlier -7.04\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHMxJREFUeJzt3X+UJWV95/H3Z2Yk2iiCM6P8GKYboiaLJivSKhjWTYDkgEHx5Jhd94yEZLPbWfC46Jq44myysieTbAhnlUQhp4MmhulVE0KM4Wz8gT+ic3bF9AgqigoxM0MHkEEElDH8mu/+UdVyp7l1f1XVrbrP/bzO6dPd99ateurWrU899TxP1VVEYGZm6VjXdAHMzKxaDnYzs8Q42M3MEuNgNzNLjIPdzCwxDnYzs8Q42K11JM1JCkkb8v//VtIFTZerSFvLJ2mPpLOaLoeNn4PdSpP0y5K+IumApLslXSXpyCFe3zOAIuKciHh/NaWtXhPlc2hbLw52K0XSW4DfA34DeCZwKjALfELSYQ2XbUOTyzdrioPdRibpCOBS4I0R8dGIeDQi9gD/hizcX59P96eSfrvjdT8taSX/+xpgK/A3kr4v6a1dlvMZSf+h4/9/L+lWSd+V9DFJsx3PhaQ3SLoNuE2Zd0q6R9IDkr4s6YUF63NILVjSOyTtzP9+qqSdkr4j6X5Jfy/pOWvLl5+97JJ0eV6+f5R0Tsc8T5D0WUnfk3SDpPesLqNLeTZJuj5f3n2SPidpXdF7Jul8SXvzMm7vsekscQ52K+PlwFOB6zofjIjvA38L/Gy/GUTE+cA+4FUR8fSIuKzX9JJeA7wd+AVgM/A54ANrJnsN8DLgJODngFcAzweOBP4t8J1+5eriArIzkuOBjcB/An5QMO3LgG8Am4DLgPdKUv7c/wa+kM/jHcD5PZb5FmCFbD2fQ7be0e09k3QScFU+v2Pz+W8ZYT0tAQ52K2MTcG9EPNblubvy56v2a8DvRsSt+XJ/B3hRZ609f/6+iPgB8CjwDODHAeWvu2uE5T5KFpbPjYjHI2J3RDxYMO3eiPjjiHgceD9wDPAcSVuBlwC/FRGPRMQu4CN9lnkMMJufDX0uim/u9Frg+oj4bEQ8DPwmcHD41bQUONitjHuBTQVt2cfkz1dtFrgib564H7gPEHBcxzR3rP4REZ8C3g28B/i2pMW8CWlY1wAfAz4o6U5Jl0l6SsG0d3cs/0D+59PJatL3dTx2SFm7+H3gduDjkr4l6W09pj2WQ9f7IUY7M7EEONitjP8HPEzWLPJDkg4HzgE+mT/0EDDTMcnRa+YzzC1G7wB+LSKO7Ph5WkT836L5RcQfRMQpwAvImmR+o2DeheXMa8yXRsRJZE1Q5wK/NES5ITuLeZakzmUcXzRxRHwvIt4SEScCrwL+i6QzV5/uMu8fzitfxsYhy2eJcLDbyCLiAbLO0z+UdLakp0iaA/6CrG34mnzSm4FXSnqWpKOBN62Z1beBEwdc7B8Bl0h6AYCkZ0r6xaKJJb1E0svy2vVDwD8DjxdMfjPwunw95smaN1bn8zOSfkLSeuBBsmaSovl0FRF7gWXgHZIOk3QaWWAXlf1cSc/N2+cfzJe3usy179m1wLmSTs9HI/0PvH9PLW94KyXv7Hw7cDlZ+NxIVqs+M2/rhSzgvwTsAT4OfGjNbH4X+G9588qv91neX5ENr/ygpAeBW8jODoocAfwx8F1gL1nzxOUF0/4m8KP5tJeSdXSuOposPB8EbgX+Dug6mqWPbcBpeTl+m+y9eLhg2ucBNwDfJzs7ujIiPpM/d8h7FhFfBd6Ql/mufB1WRiifJUD+og2z5kj6EPD1iPjvTZfF0uEau9kY5U1DP5qPRz8bOA/4cNPlsrT4yjyz8TqabNz/RrKmkgsj4qZmi2SpcVOMmVli3BRjZpaYRppiNm3aFHNzc00s2sxsYu3evfveiNjcb7pGgn1ubo7l5eUmFm1mNrEk7R1kOjfFmJklxsFuZpaYSoJd0pGSrpX09fw+2adVMV8zMxteVW3sVwAfjYjX5vepmOn3AjMzq0fpYM9vgfoK4JcBIuIR4JGy8zUzs9FU0RRzIrAf+BNJN0m6Or9t6yEkLUhalrS8f//+ChZrZmbdVBHsG4AXA1dFxMlkt0Z90hcCRMRiRMxHxPzmzX2HYZqZ2YiqCPYVYCUibsz/v5Ys6M3MrAGlgz0i7gbukPRj+UNnAl8rO18zs6btumiJlQ1zHNQ6VjbMseuipaaLNJCqRsW8EVjKR8R8C/iViuZrZtaIXRctcfJVCxxO9hW1Wx7fy1FXLbALOP3Kbc0Wro9G7u44Pz8fvqWAmbXZyoY5tjz+5Cv4V9bPsuWxPeMvECBpd0TM95vOV56amXVx7OP7hnq8TRzsZmZd3Ll+61CPt4mD3cysiz0LO3hozUX0DzHDnoUdDZVocA52M7MuTr9yGzdduMjK+lkOIlbWz3LThYut7zgFd56amU0Md56amU0pB7uZWWIc7GZmiXGwm5klxsFuZpYYB7uZWYelJZibg3Xrst9Lk3Hfr0NUdRMwM7OJt7QECwtwILvvF3v3Zv8DbGv/8PUfco3dzCy3ffsTob7qwIHs8UniYDczy+0ruL9X0eNt5WA3M8ttLbi/V9HjbeVgNzPL7dgBM4fe94uZmezxSeJgNzPLbdsGi4swOwtS9ntxcbI6TsGjYszMDrFt2+QF+VqusZuZJaayYJe0XtJNkq6vap5mZja8KmvsFwO3Vjg/MzMbQSXBLmkL8PPA1VXMz8zMRldVjf1dwFuBg0UTSFqQtCxpef/+/RUt1szM1iod7JLOBe6JiN29pouIxYiYj4j5zZs3l12smZkVqKLG/lPAqyXtAT4InCFpZwXzNTOzEZQO9oi4JCK2RMQc8DrgUxHx+tIlMzOzkXgcu5nZGIzzPu+VXnkaEZ8BPlPlPM3MJt247/PuGruZWc3GfZ93B7uZWc3GfZ93B7uZWc3GfZ93B7uZWc3GfZ93B7uZWc3GfZ9334/dzGwMxnmfd9fYzcwS42A3M0uMg92GNs4r6MxseG5jt6GM+wo6Mxuea+w2lHFfQWdmw3Ow21DGfQWdmQ3PwW5DGfcVdGY2PAe7DWXcV9BZO7kDvd0c7DaUcV9BZ+2z2oG+dy9EPNGB7nBvD0XE2Bc6Pz8fy8vLY1+umZU3N5eF+Vqzs7Bnz7hLM10k7Y6I+X7TucZuZkNxB3r7OdjNbCjuQG8/B7uZDcUd6IdqY0dy6WCXdLykT0u6VdJXJV1cRcHMBtXGHStl7kB/Qls7kkt3nko6BjgmIr4o6RnAbuA1EfG1ote489SqsvYWB5DVHqc1aGy8xt2RPLbO04i4KyK+mP/9PeBW4Liy8zUbhG9xYE1qa0dypW3skuaAk4Ebuzy3IGlZ0vL+/furXKxNsbbuWDYd2tqRXFmwS3o68JfAmyLiwbXPR8RiRMxHxPzmzZurWqxNubbuWDYd2tqRXEmwS3oKWagvRcR1VcxzErjTrnnD7FjeXla11nYkR0SpH0DAnwHvGvQ1p5xySky6nTsjZmYisr7w7GdmJnu87XbujJidjZCy35NQ5l4GWZ9J3l5mq4DlGCBjqxgVczrwOeArwMH84bdHxP8pek0Ko2KKesMhO2rv2NGCo3YX0zqKxJfBWwoGHRXje8WMaN26rN5XpK1hOa0BV7S9JDh48MmPm7WR7xVTs36dc20dcjeto0jcyXoo9zfUoy3vq4N9RN067dZqY1hOa8C1dfRCE9p6teSka9P76mAfUWdveJE2huW0BlxrRy80wBd11aNN76vb2CswaR2SS0vZh23fvuzg09aOXquH+xvqMY731W3sYzRptcFt27KO0oMHs99tLafVo+nmuLa0Q1et6fe1k4O9Ig7LeqQaAk0aV3Nct23XpnboqrWqmXOQwe5V/zRxgVJqF+VMA19UVJ+694eibbdx46GPrf7Mzla7/KbU/b4yrguURlG2jX3YNuJJawO3zLSOuU9Brwv4unH7/mCSbWPvdSpXdNrept5qG9y0jrlPwbDbqI0jyCbZxAV7UUhffHFx4LcxINx23F+bOqNSVsdnsWgbbdzYonbolA3SXlP1T5k2dql7G13Rz+xs9tOmdr1+bcfuD8iM2sbu929wdfVj9JrvsNvH2/MJDNjGPnHBXhTSRT9Ss51w3T6UvQ4009phWLTzdj6+cWP247s4VqfOSs+gd93sNY2356GSDfZuG7pfjX31deM+6hd9KHsdhNp2djEOg+y8g+7gvQ78017b66boDFiqf9mDbNNp3B96STbYI3rXett0ZC8q4/r1xR/WJne0pgyy8w66g/drqmv6M9G2ZoUmg3OQZU/j/tBL0sG+qtdOPMhp+1pV7XSDHHiKairTWEMZZOcddAcf5IDf1r6VaSvTINt0GveHXqYi2Is2+saNw39Yq/qAD9JU1NnW3q1NuW07/6BGPTCWqbF3vp+rZej3/jdV22trSDV1FjHI+zHJ+0MdpiLYq7y6raqdbpCaeoqjAMrsgKO2sRdN3++MqakgnYRmhXF+9gb9zEzi/lCXqQj2iO4bfZQdqKqdrlfzUIofyqpCdJDRL8Muq221vbbW2Fc1MQzXQx+HM9ZgB84GvgHcDryt3/R13ytmlB2o7hp7W3beKpVt9ui2kw4SxsMchNsUBG070KzV9mG4VZShTZ+HUYwt2IH1wD8AJwKHAV8CTur1mrqDfZQPQJ1t7G3aeatUpqOyTDNa1QfPKsZbV7mspvQ6YLahwlJ2KGsK++Y4g/004GMd/18CXNLrNaMEe12nbMNeAFNHWSdVmaGFgxwUimrjVe6gVY6hH4c6P1u9wrtsU2UV5S47lLUNB6eyxhnsrwWu7vj/fODdvV4zbLDXtWO1aYedRGVqUKPcGqJTVQFX5YicutX9ee01/zKhWFW5yw5lnYTO637GGey/2CXY/7DLdAvAMrC8devWoVamriNtCkfwXuo+cyizw1Y5VLWMQTq723LR0zg+r2WH4XZ7/aDl7vd5Ldunk8L+nlRTTF1H2kk4go/a/juus5FRDx69yjfOpqxBhqcWtfsPGg5VrU/Tn9deob/6Pq4tY79baHTOe5gDx6jbYdLP0McZ7BuAbwEndHSevqDXa4YN9mmtsZdp/52Eb6ppQ1/EILXAbmcRg4ZrlWFS9vNa13DFfu9Nr1tojLpuo76vbfjMlTHu4Y6vBL6Zj47Z3m/6SWtjb+rDULb9t8naXac6Oqir1K8WuHqH0FFqilVWHsrsB3XtQ4N+/qocwtq5TpMc0qNI7gKlujbiKO16bbqXRtlOyLr1q9G16VR4kBAe9vNQdfPJqPtBXWeng3z+Otvai8rd9rPntkgu2JvS5AeuTI29yk7IMgfVsiMZxqmOM7h+23Bctc662ucH6aNo8owiNQ72irT9ftVlOyHrPmMZpEY3aZ3Vg7x2tclpdf2Kts+4wqyuCkq3dVhd3zLv37Q0rQzLwV6RqttIh/3gVhHOvebdL1jKrv8k1djL6Nfk1C3s+rXZVxlsdR5EHMjj42DvYZgPYlU7RNOnmqOOL67iisNJaWMvY5QDWF3j4/sNS3QATy4He4FRAraKHaLJttaidR6kaWTUGnvRSJg2joqpwihNTnWczTRdgbB6OdgLNNUZ2qvmW/fOWLTOg4wvHvVAOG3hMkpIDzIGfNj+h1RGl/jsojsHe4GmbmbUa4ere2fsVZss6vgqU8tOJVyGMWqT06jj44s0fXVqFaaxYjAoB3uBQUOnW4DXdYFI3Ttjr3XuDJZ+B4CyI2EmKVxGUeZCrKrCLIWDagrrUBcHe4EyQwjLXqZfVNuv+4NcZvTLKOXxjjmaKpofqhx+2JRprRgMwsHew6hXwRX9lP3AjePUs986VznefJpOpdvYFtzrLGwStoMrBsUc7CU0cZl+0wFR9QiNptdnHJo+gI1aQWl7QFb9vqb0WXSwl1C0Q4z7XuHjNC3jzatUZ3D2a68fJPwmuUmjqjBu+uBbNQd7Cb0+DG05+tdRjjKdf5OmivevruAc5CDbxHfDTqLU3gMHe0ltCfBuUquFjFvbR6AM28dTdFCZts9Jt312ks9aunGwJyy1Wsi4VfX+1RWcw/bx9FqHshWUNldwOtU1kq1tHOwJS60WMm5Vvn91BN8gNfZx9PdMUo2/13s2iSODigwa7OuwibN163CP26GqfP+2bYM9e+Dgwez3tm1lSpbZsQNmZoqfn5mBK66AxUWYnQUJNm6Epz0Nzj8f5uZgaal8ObZvhwMHDn3swIHs8bbZt6/4uYjsPYLs/VpcrGY7tZmDfQJ12/FnZrLHrb+2v3/btj05tDduzP7uDKbVg8o118APfgDf+U4WYnv3wsJC+XAvCsteIdqUfgfliOy9q+rg23qDVOur/nFTTHmT0vbZVim9f+PuxG1j+3QdN1RrIwZsilE27Wgk/T7wKuARsi+y/pWIuL/f6+bn52N5eXnk5ZrZE9aty6JrLSlrIhrV0lJW8+9sjpmZaW9TxtJS1ky0d2/351dr7JNM0u6ImO83XdmmmE8AL4yInwS+CVxScn5mNqS6+lzWNgm1vX16tWlq5852N7WNQ6lgj4iPR8Rj+b+fB7aUL5IVWVrKOsbWrauug8wmX519BnV0Dtdt0g5IdSjVFHPIjKS/AT4UETsLnl8AFgC2bt16yt6i8yXratJOi228Vpsh9u3Lauo7dvhzkaJBm2L6BrukG4Cjuzy1PSL+Op9mOzAP/EIMcKRwG/vw5ua6tx2m0G5oZoMZNNg39JsgIs7qs6ALgHOBMwcJdRvNJA09M7NmlWpjl3Q28F+BV0fEgX7T2+h8UZKZDarsqJh3A88APiHpZkl/VEGZrIu2X1RjZu3Rtymml4h4blUFsd5WO8LcQWZm/ZQKdhuv1cvIzcx68b1iEuOx7maDSXlfcY09IWvHuq/eDApc0zfrlPq+4hp7QibpNquTIOUa3bRLfV9xjT0hHutendRrdNMu9X3FNfaEeKx7dVKv0U271PcVB3tCPNa9OqnX6KZd6vuKgz0hvqtddVKv0U27bt9SVfVXCzbJwZ6YSbzNahulXqOz+r9asEkOdrMufPYzPVLsT6nsfuzD8G17zawt6vpqwTqM66vxzMwmWor9KQ52M5tqKfanONjNbKql2J/iK0/NbOqldudU19jNzBLjYDczS4yD3cwsMQ52M7PEVBLskn5dUkjaVMX8zMxsdKWDXdLxwM8Cvu+dmVkLVFFjfyfwVmD89yYwM7MnKRXskl4N/FNEfGmAaRckLUta3r9/f5nFmplZD30vUJJ0A3B0l6e2A28Hfm6QBUXEIrAI2U3AhiijmZkNoW+wR8RZ3R6X9BPACcCXJAFsAb4o6aURcXelpTQzs4GN3BQTEV+JiGdHxFxEzAErwIsd6maWuqWl7JuW1q1r5zcu+V4xZmZDWFrKvmFp9cs5Vr9xCdpzv5nKLlDKa+73VjU/M7M2moRvXPKVp2ZmQ9hXcMVO0eNNcLCbmQ1hEr5xycFuZjaESfjGJQe7mdkQJuEblzwqxsxsSG3/xiXX2M3MEuNgNzNLjIPdzCwxDnYzs8Q42M3MEuNgNzNLjIPdzCwxDnYzs8Q42M3MEuNgNzNLjIPdzCwxDnYzs8Q42M3MElM62CW9UdI3JH1V0mVVFMrMzEZX6ra9kn4GOA/4yYh4WNKzqymWmZmNqmyN/ULgf0bEwwARcU/5IpmZWRllg/35wL+SdKOkv5P0kioKZWZmo+vbFCPpBuDoLk9tz19/FHAq8BLgzyWdGBHRZT4LwALA1jZ966uZWWL6BntEnFX0nKQLgevyIP+CpIPAJmB/l/ksAosA8/PzTwp+MzOrRtmmmA8DZwBIej5wGHBv2UKZmdnoyn6Z9fuA90m6BXgEuKBbM4yZmY1PqWCPiEeA11dUFjMzq4CvPDUzS4yD3cwsMQ52M7PEONjNzBLjYDczS4yD3cwsMQ52M7Ma7LpoiZUNcxzUOlY2zLHroqWxLbvsBUpmZrbGrouWOPmqBQ7nAABbHt/LUVctsAs4/cpttS/fNXYzs4rNLW7/YaivOpwDzC1uH8vyHexmZhU79vF9Qz1eNQe7mVnF7lzf/dbkRY9XzcFuZlaxPQs7eIiZQx57iBn2LOwYy/Id7GZmFTv9ym3cdOEiK+tnOYhYWT/LTRcujqXjFEBN3GV3fn4+lpeXx75cM7NJJml3RMz3m841djOzxDjYzcwS42A3M0uMg93MLDEOdjOzxJQKdkkvkvR5STdLWpb00qoKZmZmoylbY78MuDQiXgT8Vv6/mZk1qGywB3BE/vczgTtLzs/MzEoqe9veNwEfk3Q52UHi5UUTSloAFgC2bh3P/RLMzKZR3xq7pBsk3dLl5zzgQuDNEXE88GbgvUXziYjFiJiPiPnNmzdXtwZmZmPQ5BdnDKvULQUkPQAcGREhScADEXFEv9f5lgJmNknWfnEGZDf1Guf9X2B8txS4E/jX+d9nALeVnJ+ZWes0/cUZwyrbxv4fgSskbQD+mbwN3cwsJU1/ccawSgV7ROwCTqmoLGZmrXTn+q1seXxv98cbKE8/vvLUzKyPpr84Y1gOdjOzPpr+4oxh+Ys2zMwmhL9ow8xsSjnYzcwS42A3M0uMg93MLDEOdjOzxDjYzcwS42A3M0uMg93MLDGNXKAkaT/w5BsvDGYTcG+FxZkU07je07jOMJ3rPY3rDMOv92xE9P1Ci0aCvQxJy4NceZWaaVzvaVxnmM71nsZ1hvrW200xZmaJcbCbmSVmEoN9sekCNGQa13sa1xmmc72ncZ2hpvWeuDZ2MzPrbRJr7GZm1oOD3cwsMRMV7JLOlvQNSbdLelvT5amDpOMlfVrSrZK+Kuni/PFnSfqEpNvy30c1XdaqSVov6SZJ1+f/nyDpxnydPyTpsKbLWDVJR0q6VtLX821+WurbWtKb88/2LZI+IOmpKW5rSe+TdI+kWzoe67ptlfmDPNu+LOnFZZY9McEuaT3wHuAc4CTg30k6qdlS1eIx4C0R8S+AU4E35Ov5NuCTEfE84JP5/6m5GLi14//fA96Zr/N3gV9tpFT1ugL4aET8OPAvydY/2W0t6TjgPwPzEfFCYD3wOtLc1n8KnL3msaJtew7wvPxnAbiqzIInJtiBlwK3R8S3IuIR4IPAeQ2XqXIRcVdEfDH/+3tkO/pxZOv6/nyy9wOvaaaE9ZC0Bfh54Or8fwFnANfmk6S4zkcArwDeCxARj0TE/SS+rYENwNMkbQBmgLtIcFtHxGeB+9Y8XLRtzwP+LDKfB46UdMyoy56kYD8OuKPj/5X8sWRJmgNOBm4EnhMRd0EW/sCzmytZLd4FvBU4mP+/Ebg/Ih7L/09xe58I7Af+JG+CulrS4SS8rSPin4DLgX1kgf4AsJv0t/Wqom1bab5NUrCry2PJjtWU9HTgL4E3RcSDTZenTpLOBe6JiN2dD3eZNLXtvQF4MXBVRJwMPERCzS7d5G3K5wEnAMcCh5M1Q6yV2rbup9LP+yQF+wpwfMf/W4A7GypLrSQ9hSzUlyLiuvzhb6+emuW/72mqfDX4KeDVkvaQNbGdQVaDPzI/XYc0t/cKsBIRN+b/X0sW9Clv67OAf4yI/RHxKHAd8HLS39arirZtpfk2ScH+98Dz8t7zw8g6XD7ScJkql7ctvxe4NSL+V8dTHwEuyP++APjrcZetLhFxSURsiYg5su36qYjYBnwaeG0+WVLrDBARdwN3SPqx/KEzga+R8LYma4I5VdJM/llfXeekt3WHom37EeCX8tExpwIPrDbZjCQiJuYHeCXwTeAfgO1Nl6emdTyd7BTsy8DN+c8rydqcPwnclv9+VtNlrWn9fxq4Pv/7ROALwO3AXwA/0nT5aljfFwHL+fb+MHBU6tsauBT4OnALcA3wIylua+ADZP0Ij5LVyH+1aNuSNcW8J8+2r5CNGhp52b6lgJlZYiapKcbMzAbgYDczS4yD3cwsMQ52M7PEONjNzBLjYDczS4yD3cwsMf8fXa8kbimeyVYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x16bdd0185f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "std = np.std(total_data)\n",
    "mean = np.mean(total_data)\n",
    "b = 3\n",
    "outliers = []\n",
    "outlier_index = []\n",
    "lower_limit = mean - b * std\n",
    "upper_limit = mean + b * std\n",
    "\n",
    "for i in range(len(total_data)):\n",
    "    if total_data[i] < lower_limit or total_data[i] > upper_limit:\n",
    "        print('Outlier %0.2f' % (total_data[i]))\n",
    "        outliers.append(total_data[i])\n",
    "        outlier_index.append(i)\n",
    "\n",
    "plt.figure(3)\n",
    "plt.title('Outliers using std')\n",
    "plt.scatter(range(len(total_data)), total_data, c='b')\n",
    "plt.scatter(outlier_index, outliers, c='r')\n",
    "plt.savefig('B04041 04 10.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 局部异常因子(Local Outlier Factor,LOF）检测异常点\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5.0, 2.0)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = [1, 2, 12]\n",
    "np.mean(a), np.median(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "import numpy as np\n",
    "\n",
    "instances = np.matrix([[0, 0], [0, 1], [1, 0], [1, 1], [5, 0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 2\n",
    "distance = 'manhattan'\n",
    "from sklearn.metrics import pairwise_distances\n",
    "dist = pairwise_distances(instances, metric=distance)\n",
    "\n",
    "import heapq\n",
    "k_distance = defaultdict(tuple)\n",
    "for i in range(instances.shape[0]):\n",
    "    distances = dist[i].tolist()\n",
    "    ksmallest = heapq.nsmallest(k + 1, distances)[1:][k - 1]\n",
    "    ksmallest_index = distances.index(ksmallest)\n",
    "    k_distance[i] = (ksmallest, ksmallest_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "defaultdict(tuple,\n",
       "            {0: (1.0, 1), 1: (1.0, 0), 2: (1.0, 0), 3: (1.0, 1), 4: (5.0, 0)})"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k_distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All neighbourhood of point 0.00\n",
      "[0.0, 1.0, 1.0, 2.0, 5.0]\n",
      "K neighbourhood of point 0.00\n",
      "[1.0, 1.0]\n",
      "All neighbourhood of point 1.00\n",
      "[1.0, 0.0, 2.0, 1.0, 6.0]\n",
      "K neighbourhood of point 1.00\n",
      "[1.0, 1.0]\n",
      "All neighbourhood of point 2.00\n",
      "[1.0, 2.0, 0.0, 1.0, 4.0]\n",
      "K neighbourhood of point 2.00\n",
      "[1.0, 1.0]\n",
      "All neighbourhood of point 3.00\n",
      "[2.0, 1.0, 1.0, 0.0, 5.0]\n",
      "K neighbourhood of point 3.00\n",
      "[1.0, 1.0]\n",
      "All neighbourhood of point 4.00\n",
      "[5.0, 6.0, 4.0, 5.0, 0.0]\n",
      "K neighbourhood of point 4.00\n",
      "[4.0, 5.0]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "defaultdict(list,\n",
       "            {0: [(1.0, 1), (1.0, 2)],\n",
       "             1: [(1.0, 0), (1.0, 3)],\n",
       "             2: [(1.0, 0), (1.0, 3)],\n",
       "             3: [(1.0, 1), (1.0, 2)],\n",
       "             4: [(4.0, 2), (5.0, 0)]})"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def all_indices(value, inlist):\n",
    "    out_indices = []\n",
    "    idx = -1\n",
    "    while True:\n",
    "        try:\n",
    "            idx = inlist.index(value, idx + 1)\n",
    "            out_indices.append(idx)\n",
    "        except ValueError:\n",
    "            break\n",
    "    return (out_indices)\n",
    "\n",
    "\n",
    "import heapq\n",
    "k_neig = defaultdict(list)\n",
    "for i in range(instances.shape[0]):\n",
    "    distances = dist[i].tolist()\n",
    "    print('All neighbourhood of point %.2f' % (i))\n",
    "    print(distances)\n",
    "    ksmallest = heapq.nsmallest(k + 1, distances)[1:]\n",
    "    print('K neighbourhood of point %.2f' % (i))\n",
    "    print(ksmallest)\n",
    "    ksmallest_set = set(ksmallest)\n",
    "    ksmallest_idx = []\n",
    "    for x in ksmallest_set:\n",
    "        ksmallest_idx.append(all_indices(x, distances))\n",
    "    ksmallest_idx = [item for sublist in ksmallest_idx for item in sublist]\n",
    "    k_neig[i].extend(zip(ksmallest, ksmallest_idx))\n",
    "k_neig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "defaultdict(float, {0: 1.0, 1: 1.0, 2: 1.0, 3: 1.0, 4: 0.2222222222222222})"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "local_reach_density = defaultdict(float)\n",
    "for i in range(instances.shape[0]):\n",
    "    no_neighbours = len(k_neig[i])\n",
    "    denom_sum = 0\n",
    "    for neigh in k_neig[i]:\n",
    "        denom_sum += max(neigh[0], k_distance[neigh[1]][0])\n",
    "    local_reach_density[i] = no_neighbours / (1.0 * denom_sum)\n",
    "\n",
    "local_reach_density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "lof_list = []\n",
    "for i in range(instances.shape[0]):\n",
    "    lrd_sum = 0\n",
    "    rdist_sum = 0\n",
    "    for neigh in k_neig[i]:\n",
    "        lrd_sum += local_reach_density[neigh[1]]\n",
    "        rdist_sum += max(neigh[0], k_distance[neigh[1]][0])\n",
    "    lof_list.append((i, lrd_sum * rdist_sum))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, 4.0), (1, 4.0), (2, 4.0), (3, 4.0), (4, 18.0)]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lof_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
